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It is demonstrated that the signatures of the Hubbard Model in the strongly interacting regime can be 
simulated by modifying the screening in the limit of zero wavevector in Projector-Augmented Wave GW 
calculations for systems without significant nesting. This modification, when applied to the Mott insulator 
CuO, results in the opening of the Mott gap by the splitting of states at the Fermi level into upper and lower 
Hubbard bands, and exhibits a giant transfer of spectral weight upon electron doping. The method is also 
employed to clearly illustrate that the Mi and M 2 forms of vanadium dioxide are fundamentally different types 
of insulator. Standard GW calculations are sufficient to open a gap in Mi VO 2 , which arise from the Peierls 
pairings filling the valence band, creating homopolar bonds. The valence band wavefunctions are stabilized 
with respect to the conduction band, reducing polarizability and pushing the conduction band eigenvalues to 
higher energy. The M 2 structure however opens a gap from strong on-site interactions; it is a Mott insulator. 
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I. INTRODUCTION 

Strongly correlated electrons in transition metal oxides 
generate phenomena such as high temperature supercon¬ 
ductivity, colossal magnetoresistance and metal-insulator 
transitions, which offer enormous potential for new gen¬ 
erations of devices The development of density func¬ 
tional theory (DFT)^ and its applications to weakly cor¬ 
related materials such as the p-block semiconductors has 
significantly facilitated material design. Strongly corre¬ 
lated materials however, have not received such benefits 
as accurate approximations to the quantum many-body 
problem have proven elusive 

The paradigm for the description of strong electron 
correlations is the Hubbard Hamiltonian^ (equation 1), 
which exhibits a competition between hopping, given by 
the t term, and repulsion, given by the U term, which 
lies at the heart of strongly correlated systems. If the 
orbitals between which the electrons hop are localized, 
such as transition metal d- or f-orbitals, then hopping to 
an already occupied site incurs an energy penalty, U. 

H = —t + U (^) 

Sj) * 

Therefore, at half-filling a tendency towards sites being 
singly occupied, and thus insulating behavior is observed. 
This localization is a many-body effect: the U term is in¬ 
curred by electrons encountering each other, and there¬ 
fore the localization of each electron requires knowledge 
of where the other electrons are. Since the discovery of 
high temperature superconductivity, attempts to develop 
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ah initio tools to describe this have taken on increased ur¬ 
gency. To-date only two approaches have emerged which 
have been able to generate significant insight. 

Dynamical Mean Field Theory^ in its DFT-j-DMFT 
form takes input wavefunctions from Density Func¬ 
tional Theory for the bands of interest, then applies 
sophisticated Monte Garlo approaches to the interac¬ 
tions of these electrons via the model Hubbard Hamil¬ 
tonian. DMFT has provided significant insight into 
the nature of some strongly correlated systems such as 
V 203 ^ LaOi_a^Fa,FeAs, and FeSe;^ The technique re¬ 
volves around isolation of the bands of interest, and 
projecting them into real space, usually onto Wannier 
functions^ which are a natural basis for the Hubbard 
Model (indeed, this was the basis chosen by Hubbard 
himself). Gontributions to the electron self-energy from 
screening and scattering are evaluated using a Mean Field 
Approximation based on the Anderson Impurity Model, 
or cluster variationsDespite such successes, DMFT and 
its variations require considerable power and sophistica¬ 
tion when applied to real systems. Other complementary 
techniques therefore become attractive when the number 
of correlated bands in the system is large. 

The GW approximationiSiii is a many-body perturba¬ 
tion theory approach which in the last decade has been 
integrated with density functional methods to provide ac¬ 
curate ab initio calculations of the electronic structures of 
materials such as Si, GaAs etc j^^d^ This approach takes 
the bare Hartree-Fock interaction, well known to over¬ 
state the interactions between electrons, and screens it 
with the dielectric matrix, usually calculated in the Ran¬ 
dom Phase Approximation (RPA). This approach has ex¬ 
hibited significant improvements over DFT, and is now a 
standard component of most ab initio packages. How¬ 
ever, to maintain computational tract ability, the self¬ 
energy is constructed using non-interacting Green func- 
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tions and completely neglects electron scattering dia¬ 
grams, and thus significant modification is required to 
apply the GW Approximation to Mott insulating sys¬ 
tems. 

By far the most common method of “correcting” the 
RPA-based GW method is to modify the input wavefunc- 
tions and eigenvalues by calculating them using either the 
DFT+U approach of Anisimov et al^ or with hybrid 
functionalsJiang et al. employed the LDA-l-U method 
to the electronic structures of oxides of the lanthanides^ 
and the first row transition metal oxides MnO, FeO, GoO 
and NiO^ii Gonsiderable improvement was found in the 
band gaps in both studies. 

Rodl et compared self-consistent GW calculations 
using input wavefunctions from DFT+U and hybrid func¬ 
tional calculations utilizing different iteration schemes for 
the calculation of the photoemission spectrum of GuO. 
Of these, the authors found that self-consistency in both 
eigenvalues and wavefunctions vastly overestimated the 
band gap if the static screening was first approximated 
by the HSE06 functional^ or PBEd^-hU calculations. 
The authors found that an approach in which a screened 
Goulomb interaction in good agreement with experiment 
was used and held fixed achieved the best approximation 
to experimental data. 

Lanyas introduced an arbitrary (attractive) on-site ad¬ 
dition to the local potential to GW calculations of 3d 
metal oxides (including GuO) in an attempt to both ob¬ 
tain band gaps with better agreement with experiment, 
and fix the incorrect band ordering generated by using 
hybrid functional input to GW calculations. Again con¬ 
siderable improvement was found for band gaps and band 
ordering. We discuss the basis of the DFT+U and hybrid 
functional input approach in comparison to the method 
developed in this work in Section IIG. 

Gatti and Guzzo^ applied the “GW+G” method (the 
G stands for Gumulant), which is obtained from a de¬ 
coupling of the elements of the exact one-electron Green 
function in the Dyson equation^ to the study the effect 
of satellites on the electronic structure of SrVOs. This 
approach successfully renormalized the V 3d bands and 
satellites near the Fermi level from first principles. This 
renormalization gave good agreement with experiment, 
without the use of model Hamiltonians. 

Recently, some attempts to apply GW to strongly 
correlated systems have focused on combining it with 
DMFT, however the significant theoretical and com¬ 
putational complexity has so far limited this approach 
to model Hamiltonians adatoms on surfaces^ and 
SrVOs 

In this work, we take a different approach. Rather 
than adjust the input wavefunctions, we approximate the 
effect of on-site repulsion by partially unscreening the 
Goulomb interaction in the limit of low wavevector in the 
GW calculation. This approach simulates the scattering 
resulting from the Hubbard U term by mimicking the 
closing of the polarization bubbles by electrons on other 
sites. We refer to this technique as “Partially Screened 


GW,” or PS-GW. 


II. METHODS 
A. Assumptions 

The approach taken in this work is based around the 
following assumptions: i) the GW Approximation based 
on the screened interaction calculated in the RPA gets 
the electronic structure “mostly” correct. This assump¬ 
tion infers that RPA-based GW calculations can be “cor¬ 
rected” for the effects of scattering in the self-energy, ii) 
the most significant contributions to the self-energy of 
electrons near the Fermi level (E^?), come from inter-site 
hopping in the systems studied in this work, iii) these 
transitions manifest in the low q limit of the dielectric 
response, and thus the transitions which create double 
occupancies manifest in this limit, iv) correlations will 
heavily suppress transitions for high frequencies, such 
that the static limit of X^(q, cj) is a more accurate repre¬ 
sentation of the response of real systems. 

Of these three assumptions, ii) and iv) are the most 
difficult to justify. Assumption iii) follows from ii) and 
is supported by calculations; if ii) holds then calculating 
e(q, cj) on different q-point grids should give the same 
results, and this is indeed the case. For all systems stud¬ 
ied in this work, the dielectric response only depends on 
the low q transitions. Ghanging the number of q points 
leaves the response invariant. Justification of assump¬ 
tion i) follows from the data itself. If RPA-GW can be 
corrected (using sound theoretical arguments) such that 
agreement with experimental results is achieved, then the 
approach is justified. We take the same approach to the 
justification of assumption iv), which has the added ben¬ 
efit of significantly reducing the computational resources 
required. 

Assumption ii) depends to a large extent on the 
method used to generate the input wavefunctions. The 
two Mott systems studied here, GuO and M 2 VO 2 , are 
both Monoclinic. This Monoclinic structure results from 
the adoption of a charge density wave which stabilizes 
the oxygen states, reducing their energies with respect to 
the metal d-states at the Fermi level. This stabilization 
arises from a change in nuclear potential, and is thus well 
reproduced by DFT. For these systems, DFT provides an 
adequate starting point for the method described here. 
For other systems, such as Gu 2 0^ this is not the case, 
and the method used to generate the input wavefunctions 
must be modified accordingly, such as by employing the 
DFT+U method^ or Hybrid Functionals 

Thus the overall premise of this work is that: in ab 
initio calculations of Mott systems, the creation of double 
occupancies results from the low q limit of the dielectric 
matrix, and that modification of the RPA screening in 
this limit can reproduce the signatures of Mott systems. 
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FIG. 1. Feynman diagrams of relevant interactions in the systems studied, a) Polarization bubble of the type summed over to 
generate the screened interaction in the Random Phase Approximation, b) scattering vertex F(k, k') resulting from the inter¬ 
action of the polarization bubble with another momentum state, c) as per b) but with a large overlap shifting the characteristic 
frequency of the bubble to higher cj, d) scattering vertex with Hartree interaction at later time, e) low q limit of the vertex 
and f) low q limit, and limit as oj ^ oo of the vertex. 


B. GW Calculations 


ciprocal space, (x^(cc;)) as per)^ 


The GW approximation is encapsulated by the Hedin 
equations^ which constitute a self-consistent approach 
to evaluating the electron self-energy: 


^0 = -iGGT 

(2) 

(5E 



(3) 

W = e-^iy 

(4) 

O 

1 

1—1 

II 

(5) 

S = iGWT 

(6) 


X°(q,w) = ^ 'p+q /np) 

^ GG' nn' 

^ T ^n'p-|-q ^np T ^n'p-|-q] 


In the VASP implementation of the GW approximation 
used in this study, the polarizability is combined with the 
Coulomb interaction into the dielectric matrix (equation 
(5)) as per>22 


where G is a single particle Green function, is the 
irreducible polarizability, F is a Vertex Function, W is 
the screened interaction, e is the quantum dielectric ma¬ 
trix, and X is the self-energy. In ab initio calculations, 
X is used in place of the DFT exchange-correlation en¬ 
ergy in a “quasiparticle” Hamiltonian to generate band 
eigenvalues ^22 The polarizability matrix, x = —iGGF, 
is usually evaluated by setting F=1 (unless excitonic ef¬ 
fects are important, in which case F is approximated by 
solving the Bethe-Salpeter equation^), and the Green 
functions used are non-interacting i^^i^^i^^ Standard per¬ 
turbation theory provides a computational form of this 
independent-particle polarizability matrix weighted by 
the Hartree-Fock overlap integrals for the vertex in re- 


Cq(Gr, G^, Cj) — fc,G' “ 


47re^ 


|q + G||q + G^ 


XXq(G, G',Cj) 


( 8 ) 


This dielectric matrix screens the bare Coulomb in¬ 
teraction to generate the screened interaction (equation 
(4))^ 


Wq(G,G',w) = 47re^ ' 


|q + G 


■£q (G, G',w) 


q+G'l 


(9) 

The computational version of equation (6) is then given 
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FIG. 2. a) Hartree-Fock calculation of CuO b) Unmodified frequency-dependent spin-resolved Go Wo calculation of the states 
near Ef of CuO, c) static PS-G3W3 calculation of CuO, d) static PS-G3W3 calculation of CuO doped with one electron, e) 
DFT (gray lines) and static PS-G3W3 (black filled circles) band structures of CuO, f) DFT (gray lines) and static PS-G3W3 
(black hlled circles) band structures of electron doped CuO. PS-GW derived eigenvalues are fitted with blue splines as a guide 
to the eye. 


by>^ 




nknk — 
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_1_> 

ui-Lo' - e„/k-q + i7?sgn[e„/k-q - l4 ' 


( 10 ) 


Thus the two vertices of the screened interaction are 
split, with the photon propagator terms appearing in 
equations (8) and (9), while the Hartree-Fock overlap 
integrals are the angle-bracket terms contained in equa¬ 
tions (7) and (10). In the VASP implementation of 
the GW approximation the bare Coulomb kernel is sub¬ 
tracted from the screened interaction to make the fre¬ 
quency integral in equation (10) well-behaved. This term 
is then added back into the self-energy 

irq(G, G', co') = Wq(G, G', w') - G') (11) 


E(cij)72k,nk — E(C(j)72k,nk T (V^nk |iV^nk) 


From equations (4) and (5) we see that the screened 
interaction W (q) is bounded from above by the Hartree- 
Fock interaction, and thus in real systems ranges from 
small values, when the dielectric response e(q) is large, 
up to this Hartree-Fock limit. 


C. Screening and Scattering 


Figure 1 lists some of the relevant Feynman dia¬ 
grams for the processes under consideration in this 
study. Figure la illustrates a polarization bubble, which 
are summed over in the Random Phase Approximation 
(RPA) to generate the RPA self-energy (equations 7-10). 
It represents an electron in momentum state k emitting 
a photon of wavevector q at time leaving it in momen¬ 
tum state k — q. This photon is absorbed by state p, 
promoting it to state p + q, leaving a hole in state p. 
At time 1/ state p + q emits a photon of wavevector q, 
decaying back into state p, which is absorbed by state 
k — q returning it to the original momentum state k. 

The RPA polarizability x(q, cj) gives the amplitude for 
bubbles of this type to occur (this is equation (7) without 
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the overlap integrals): 




E 

pnn' 


/n^p+q 
^ ^n'p+q 


fnp _ 

^np E 


(13) 


This will go “on-shell” when uj — e^'p+q + Cnp ~ 0. For 
a metal with many states available near the Fermi level, 
Cn'p+q — ^np ~ 0 for many transitions, and so uj is small. 
This corresponds to t' — t ^ oo, and therefore the polar¬ 
ization bubbles are long-lived, and the electrons are well 
screened. Taking this RPA interaction as a starting point 
(Assumption i), we can explore the effect of strong corre¬ 
lations by adding an interaction with another momentum 
state, k'. 

If the transitions between the filled and the low-energy 
available states correspond to inter-site hopping (As¬ 
sumption ii)), then the polarization bubble results in a 
double occupancy. If the Hubbard effective U (equation 
(I)) is large with respect to the t term, then the am¬ 
plitude for another momentum state to interact with the 
bubble, and close it will be large. The first order approxi¬ 
mation to this scattering process^ r(k, k'), is illustrated 
in Figure Ib (exchange variants, in which each interac¬ 
tion z/(q) results in exchange are also possible but not 
pictured). The second interaction, z^'(q), is given by: 


^'(q) 


/ #r / ci3r'</>;(r)0^,_^q(r')0p+q(r)0k' (r') 
|r-r'| 


(14) 


For interactions between momentum states con¬ 
structed from d-orbitals, |r — r'| is small, and thus the 
interaction is strong. Therefore, this large amplitude will 
close the polarization bubble at some time if' < if. The 
frequency at which the scattering vertex goes “on-shell” 
shifts to higher cj, which is illustrated in Figure Ic: the 
large amplitude for the interaction between p + q and 
k' reduces the characteristic time the bubble stays open, 
shifting t" closer to t. 

This shift to higher frequency results in a stronger in¬ 
teraction between states k and p at longer times. Figure 
Id illustrates this using a Hartree bubble (although an 
exchange interaction is also possible). Since the scatter¬ 
ing process closes the bubble at time t", a bare Coulomb 
interaction is now possible at time t'", whereas without 
scattering the polarization bubble would still be open, 
and this interaction would be screened. Since the proba¬ 
bility of the scattering vertex occurring is large due to the 
on-site overlap, diagrams such as Figure Id may domi¬ 
nate those of Figure la, resulting in more interactions 
between the momentum states. 

Computing the RPA dielectric function for CuO and 
M 2 VO 2 using a single q point, and comparing this to 
a calculation using the full q grid reveals that the di¬ 
electric response is dominated by the low q transitions 
(note that this does not imply that the self-energy has a 
similar dependence), due to the photon propagator. The 
results are identical to the precision used in these calcu¬ 
lations. Therefore, if the method used to generate the 
input wavefunctions is correct, but does not account for 


the on-site interaction (e.g. DFT), these low q transi¬ 
tions correspond to inter-site hopping, generating double 
occupancies. Since the dielectric response takes the form 
of a polarization bubble multiplied by a photon propaga¬ 
tor, the creation of polarization bubbles which generate 
an energy penalty in the form of the Hubbard U term will 
be found in the low q limit, and the scattering diagram 
which closes the bubbles will be the low q limit of Fig¬ 
ure Id, which is pictured in Figure le. From momentum 
conservation, if the photon propagator gives a ^ 
dependency of the dielectric response, then the scatter¬ 
ing vertex, r(k, k'), which “corrects” the RPA screening. 
Figure Ib, must have the same q dependence at first or¬ 
der as the dielectric response. We make the ansatz that 
in the absence of significant Fermi surface nesting, the 
scattering which modifies the RPA and leads to Hub¬ 
bard physics manifests in the low q limit. In addition, 
this suggests that in the limit of q 0 the momentum 
states entering the scattering vertex are unchanged, i.e. 
forward scattering. 

From the preceding argument, the consequences of on¬ 
site interactions from Mott systems will manifest as a 
change in the frequency dependence of the RPA screen¬ 
ing. As mentioned in section I, the most common method 
of correcting this discrepancy is use the DFT+U ap¬ 
proach, or a hybrid functional such as HSE06 which 
mixes in exact exchange to calculate the input wavefunc¬ 
tions. These shift the conduction band eigenvalues to 
higher energy, thus increasing the frequency of the polar¬ 
ization response. That is, these modify e^'k+q, such that 
Cn'p+q — > 0 (cquatiou 13), increasing the frequency 

at which the process goes on-shell, and un-screening the 
low frequency interactions. However, there is another 
possibility. Rather than adjusting the eigenvalues, the 
frequency at which a bubble goes on-shell can be ad¬ 
justed. 

Strong correlations in which electron localization is ob¬ 
served experimentally suggest that there is a considerable 
shift of the frequency dependence of the bubbles. Thus 
if the Hubbard U term is large, corresponding to large 
amplitudes for the interaction at t", the characteristic 
time of the polarization bubble will approach zero. This 
means that the Hartree-Fock interaction is effectively un¬ 
screened. A schematic of this is presented in Figure If, 
where the instantaneous closing of the polarization bub¬ 
ble is effectively a Hartree interaction (again, exchange is 
also a possibility), followed by further unscreened inter¬ 
actions. Thus the scattering vertex, (Figure Ib) which 
is a function of two single-particle Green functions, can 
be represented by one Green function. This reduces the 
computational load to that of a standard GW calcula¬ 
tion. 

Therefore, if on-site interactions are generated by the 
low q limit of e(q, cj), then replacing the RPA response 
with the bare Hartree-Fock interaction for on-site inter¬ 
actions will simulate strong correlations. How this is 
achieved is detailed below. 
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D. On-site Interactions in the Projector Augmented Wave 
Method 

In the PAW method^ utilized in the Vienna Ab Ini¬ 
tio Simulation Package (VASP) code^ the all electron 
wavefunctions IV^nk) are expanded as per>^ 

It^nk) = It^nk) V ^ ~ |0i)) (Pi iV^n/c) (1^) 

i 

The IV^nk) are the pseudowavefunctions, related to the 
cell periodic part of the wavefunctions through: 

I^„k)=e*'^’'|«nk) (16) 

and the are expanded in plane waves. The partial 
waves 1), are solutions of the radial Schroedinger equa¬ 
tion for a reference atom, while the pseudo-partial waves 
(|0i)) are equivalent to the \(j)i) outside a core radius Tc, 
and the projector functions are dual to the partial waves, 
= ^ij‘ The subscript i thus denotes the atomic 
position Ri, angular momentum k and and the band 
index. Equations (7) and (10) require matrix elements of 
the form: 

(V’„'k+q|e-^'»'-|^nk) (17) 

to be calculated. Following Gajdos et alM this can be 

written: 

y^B„/k+q,nk( r)d^r (18) 

where: 

^n'k+q,nk(r) =6^ '0^/]^^q(r)'0nk(l*) (19) 

Inserting the PAW expansion (equation (15)) this 

becomes)^ 

B„/k+q.nk(r) = «;,k+q(r)«„k(r)+ 

^ ^ ('Wn'k+q |Pik+q) ('iink+q |Pik+q) 

X [(pi {r)(pj (r) - {r)^j (r)] (20) 

In the low q limit the exponential is expanded as peri^ 

giq(r-R0 _ ^ ^ ^21) 

The cell-periodic parts of the wavefunctions in k-space 
are also expanded to first order around the valence wave- 
functions: 

"^nk+q ~ "^nk T Q^k'^nk T ^(Q ) 

The low q limit of the exchange charge density then 
becomes)^ 

l™('^n'k+qk“'‘^n'0nk) = |q| (q/^n'kl^nk) (23) 


where q = q/|q| and^ 

|/^nk) = (1 T ^ ^ |Pik)Qijf (Pjk I) I Vk'^nk)“i“ 

i 

V \Pi\^)Qij r^i)|'^nk)) '^( ^ ^ |Pik)'^ijf (Pjfk I'^nk) ) 

ij ij 

(24) 

with 

Qij = [ [(t)i{r)(l)j{r) - ^i{T)^j{T)]d^r (25) 

Jq(paw) 


Tij= / (r-R^)[0i(r)0j(r)-0i(r)0j(r)]c^^r (26) 

Jn(PAW) 

By setting |Vk?^nk) = 0 in equation (23), the 

planewave terms are eliminated, leaving only the aug¬ 
mentation sphere contributions. All of the cell pe¬ 
riodic functions occur in overlap integrals with the 
projector functions, which are just the expansion co¬ 
efficients for the all electron wavefunction (equation 
15). The (2n'k+q|?^nk) drop out due to orthogonality, 
and the only other cell-periodic terms are of the form 
(Vk?^n'k+q|---|?^nk), wffich are set to zero. To see this 
more clearly, we can write the charge density at a point 
r the for two orbitals a and b as^l 

V’a(r)V’6(r) = riabir) = riabir) - n^ft(r) + nl^ir) (27) 

where na6(r) comes from the planewave expansion in 
equation (16), and the other two terms with the super¬ 
script 1, n^ 5 (r) and n^^(r) are one-center terms coming 
from the \(j)i) and \(j)i) respectively. They are only eval¬ 
uated on the radial PAW grid^ and thus only overlap 
when they correspond to the same atomic site, R^ = Rj. 
Setting |Vk?^nk) = 0 in equation (23) is equivalent to 
eliminating the first term on the right-hand side of equa¬ 
tion (27), leaving only the one-center terms. This reduc¬ 
tion means that the interaction is non-zero only between 
the atomic-like wavefunctions inside the augmentation 
spheres on the same site, i.e. the interaction is on-site. 

In the implementation of the GW approximation used 
in this study^ the dielectric matrix of equation (8) is set 
equal to the unit diagonal whenever terms inside atomic 
spheres are evaluated, resulting in a bare Hartree-Fock 
interaction. Thus, by eliminating the plane wave terms 
in the low q limit, a very strong (Hartree-Fock) on-site 
only interaction replaces the screened interaction in the 
self-energy. Unscreening the interaction in this manner 
will therefore simulate the behavior of the scattering ver¬ 
tex of Figure lb, as the polarization bubbles are largely, 
although not completely, eliminated as is explained be¬ 
low. 

From equations (18) and (19) it is also clear that the 
magnitudes of the bare Goulomb terms are controlled 
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in part by the derivatives of the cell periodic parts of 
the wavefunctions, arising from the first order expansion 
(17). These derivatives are evaluated in gapless systems 
using second order perturbation theory, which creates is¬ 
sues for self-consistency, as once the first round of quasi¬ 
particle shifts are evaluated the DFT Hamiltonian is no 
longer valid. Such gapless systems also include Mott in¬ 
sulators if DFT is used to generate the input wavefunc¬ 
tions. It is obvious however, that setting these terms 
to zero, |Vk?^nk) = 0 , will both remove this obstacle to 
self-consistency, and reduce the magnitudes of the over¬ 
lap integrals which weight the interaction (' 0 nk|^cc|' 0 nk)- 
Therefore, the removal of the | Vk^nk) terms from the cal¬ 
culation has two effects. They both completely unscreen 
the low q interactions, by setting the dielectric matrix 
to the unit diagonal, and also reduce the magnitudes of 
the bare Coulomb terms added into the self-energy. This 
effectively replaces the screened interaction at low q with 
a stronger interaction, however in which small amount of 
screening is still present. Note that when the self-energy 
is evaluated as per equation ( 10 ), the full PAW wavefunc¬ 
tions are used, which are orthonormal. 

The limitation of this approach to long range interac¬ 
tions in this work renders the method applicable only to 
low symmetry Mott systems which do not exhibit sig¬ 
nificant Fermi surface nesting at finite wavevector. If 
nesting is present, then the polarizability x(q) will ex¬ 
hibit a peak at the nesting wavevector, corresponding to 
a significant amplitude for particle-hole formation. Since 
our method does not penalize pair bubbles at finite q, 
the interactions will be underestimated by the RPA cor¬ 
relations, as usual. However, despite this, the results for 
low symmetry systems are in general quite illuminating. 


E. Computational Details 

The computational scheme used for CuO was as 
follows. The experimentally determined structural 
parameters^ were used as input to density functional 
theory^i^ calculations on 8 x 8 x 6 Monkhorst-Pack k- 
space grids, using the Generalized Gradient Approxima¬ 
tion approach to exchange and correlation of Perdew et 
al. (PBE)^ and the Brillouin zone integration approach 
of Bloechl et No initial spin ordering was assumed 
in all calculations. GW calculations were performed us¬ 
ing the implementation of Shishkin and Kress o^^i^^ in 
the Vienna Ab Initio Simulation Package (VASP)^ in 
either fully frequency dependent, or static {uj = 0 ) modes 
using 256 bands. The frequency dependent calculations 
were performed as one-shot GqWq calculations, while the 
Partially-Screened GW calculations utilized three self- 
consistency iterations (G 3 IF 3 ). An energy cutoff of 200 
eV was used for all GW calculations. 

The Ml and M 2 structures used were those of 
Anderssor4S and Marezio et al^ respectively. The Mi 
structure was first relaxed to a ground state using GGA 
DFT (PBE)J^ The M 2 structure was not relaxed, as 


DET underestimates the correlation energy^ and thus 
reaches an incorrect ground state. 6 x 6 x 6 and 4x6x6 
Monkhorst-Pack k-space grids were used for the Mi and 
M 2 structures respectively and the GW calculations were 
performed again with VASP using 256 bands, after first 
calculating input wavefunctions using DET with PBE 
GGA functionals. The static PS-GW calculations of the 
M 2 structure utilized four self-consistency steps (G 4 W 4 ) 
and 256 bands. 

Gonvergence tests for both CuO and M 2 VO 2 are 
presented in the Appendix. Eor both systems conver¬ 
gence was achieved in a relatively small number of self- 
consistency steps; 3 for CuO and 4 for M 2 VO 2 . 


III. RESULTS AND DISCUSSION 
A. Application to CuO 

Eigure 2 details the application of this technique to 
the Mott insulator CuO 4^ Eigure 2 a illustrates the den¬ 
sity of states (DOS) near the Eermi level of a Hartree- 
Eock calculation, and as is commonly observed, while 
the Hartree-Eock approach does open a gap, it is over¬ 
estimated, and predicts the ground state to be Eerro- 
magnetic, again contradicting experimental data. Eig¬ 
ure 2 b presents a standard (unmodified) spin-resolved 
Go Wo calculation of the states at the Eermi level. As is 
expected from the independent particle-RPA approach, 
the non-interacting Green functions and neglect of scat¬ 
tering vertices over-screens the Hartree-Eock interaction, 
resulting in metallic behavior. In fact, very little dif¬ 
ference is exhibited between the Go Wo and DET calcu¬ 
lations (a comparison of the DET and standard Go Wo 
band structures is presented in the Supporting Informa¬ 
tion). PS-G 3 W 3 calculations however reveal a different 
story (Eigure 2 c). Clear splitting of the states at the 
Eermi level into two characteristic peaks, the upper and 
lower Hubbard bands (UHB and LHB labels on Eigure 
2 c) is exhibited (although the use of the static limit does 
broaden the lower lying oxygen bands, see Supporting In¬ 
formation), and an excitation gap of approximately 1.1 
eV has opened. 

The magnitude of this gap compares relatively favor¬ 
ably with the experimentally determined values of 1.4- 
1.7 eV^ and the overall shape of the DOS is in reason¬ 
able agreement with Rodl et apart from the afore¬ 
mentioned broadening. However there is another exper¬ 
imental signature of Mott systems which must be sim¬ 
ulated in order for the technique to be considered an 
accurate reproduction. When doped with electrons or 
holes, Mott systems exhibit a giant transfer of spectral 
weight which clearly illustrates the failure of band theory 
for these systems 4^ In a system well-described by band 
theory, if an electron is doped into the conduction band, 
the Eermi level shifts up, and the conduction band inter¬ 
sects the Eermi level with minimal change in dispersion. 
In a Mott system, adding a small number of carriers is 
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FIG. 3. a) Left panel: Mi VO 2 band structures calculated using DFT (gray lines) and unmodified frequency-dependent GqWq 
(black filled circles fitted with blue splines), right panel: corresponding DFT (gray line) and GqWq (blue filled curve) densities 
of states, b) charge density isosurface of the valence band of Mi VO 2 , and c) charge density isosurface of the conduction band 
of Ml VO 2 . Both isosurfaces perspectives correspond to the (Oil) plane. 


not expected to significantly affect the t/U balance, and 
the system is still expected to be gapped. Therefore, 
any previously empty state in the upper Hubbard band 
which is filled upon doping must then cross the gap to sit 
in the lower Hubbard band, which significantly changes 
band dispersion. This effect was clearly observed in re¬ 
cent photoemission experiments on TiOCl4^ Figure 2 d 
presents static PS-G 3 LF 3 calculations of CuO doped with 
one electron. When compared to Figure 2c, it is clear 
that spectral weight has shifted from the upper Hubbard 
band to sit at the leading edge of the lower Hubbard 
band. Figures 2e-f illustrate how this occurs. Figure 2e 
is a comparison of the DFT and static PS-G 3 IF 3 band 
structure calculations. We use DFT for this comparison 
rather than standard GW data as the standard GW data 
is virtually identical to the DFT data (see Supporting 
Information), and the lower computational cost of DFT 
allows us to use much higher k-space resolution. Figure 
2 f presents the same data for the electron doped struc¬ 
ture. The calculations of the undoped structure reveal 
that the effect of unscreening the low q interactions is 
to split the spectral weight at the Fermi level. As noted 
in numerous DMFT calculations, band crossings at the 
Fermi level still existElectron doping (Figure 2f) 
fills the lowest lying states in the upper Hubbard band, 


which Figure 2e indicates are around the F point, drop¬ 
ping them onto the leading edge of the lower Hubbard 
band, at approximately -1 eV. 


B. Application to Vanadium Dioxide 

Turning now to the technologically significant prob¬ 
lem of the natures of the insulating phases of vana¬ 
dium dioxide^-—, Figure 3 presents GqWq calculations 
on the Ml form. This structure undergoes an insulator- 
metal transition at ^ 340 K as it spontaneously changes 
from the monoclinic P2i/c structure to the tetragonal 
P 42 /mnm form^^ Figure 3a presents a comparison of 
the DFT and GqWq bands, and the respective densities 
of states and the data clearly illustrates that the GqWq 
calculations result in splitting of the bands with respect 
to the DFT calculation, with the empty conduction band 
simply shifting upwards, with minimal change in disper¬ 
sion. The gap magnitude of ^ 0.7 eV is in excellent agree¬ 
ment with the experimental value, which is also ^0.70 
eV— and the scCOHSEX-GoWq calculations of Gatti et 
al^ 

Eigures 3b-c present charge density isosurfaces of the 
valence (Eig. 3a) and conduction (Eig. 3c) bands in the 
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FIG. 4 . a) Left Panel: M2 VO2 band structures calculated using DFT (gray lines) and unmodified frequency-dependent Go Wo 
(black filled circles, fitted with blue splines), right panel: corresponding DFT (gray line) and Go Wo (blue filled curve) densities 
of states, b) charge density isosurface of the valence band of M2 VO2, c) charge density isosurface of the quasiparticle peak of 
M2 VO2, d) charge density isosurface of the conduction band of M2 VO2, e) projected densities of states of the Peierls chain 
vanadium atoms of M2 VO2, f) projected densities of states of the antiferroelectric chain vanadium atoms of M2 VO2. 


(1,1,0) plane (the “conduction band” is the first peak 
above the Fermi level in the Density of States). As ex¬ 
pected from the well-known Peierls distortion of the Mi 
structure, the pairing of the vanadium nuclei results in 
bonding density between the nuclei, while the conduction 
band consists of the corresponding anti-bonding states, 
thus confirming that the gap in Mi V02 opens via bond- 
ing/antibonding splitting. The magnitude is significantly 


underestimated by DFT however. Figures 3b-c indicate 
that Peierls pairing produces an increase in the inter¬ 
vanadium local potential, which stabilizes bonding wave- 
functions with respect to conduction states, shifting the 
conduction band eigenvalues to higher energy. 

Figure 4 illustrates that this is not the case for the M 2 
form. M 2 vanadium dioxide also undergoes an insulator- 
metal transition, although at slightly higher temperature 
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FIG. 5. a) Left Panel: M2 VO2 band structures calculated using DFT (gray lines) and static PS-G4W4 (black circles, fitted 
with blue splines), right panel: corresponding DFT (gray line) and static PS-G4W4 (blue filled curve) densities of states, b) 
static PS-G4W4 projected densities of states of the Peierls chain vanadium atoms of M2 VO2 and c) static PS-G4W4 projected 
densities of states of the antiferroelectric chain vanadium atoms of M2 VO2. 


(353 K— ), coincident with a structural transition from 
monoclinic C2/m to the same tetragonal P42/mnm struc¬ 
ture as the Ml form;^ However, the monoclinic form 
differs significantly in structure. In the Mi form all of 
the vanadium atoms form Peierls paired chains running 
down the monoclinic a-axis, which experience a slight an¬ 
tiferroelectric twist that has components in both the b- 
and c-axes. The M 2 form however, has two distinct chain 
structures (a comparison of this with the Mi structure 
is presented in the Supporting Information). One half of 
the vanadium atoms form a Peierls paired chain, how¬ 
ever this chain is not antiferroelectrically distorted, but 
rather the vanadium atoms are collinear. The remaining 
vanadium atoms form an antiferroelectrically distorted 
chain, however one in which the inter-vanadium spacing 
is uniform {i.e. no Peierls pairing). This structure has 
been regarded as a Mott insulator since the 1970s, due to 
experiments by Pouget and co- workers^ who used 
NMR Knight shifts to resolve the two vanadium environ¬ 
ments in Cr-doped VO 2 . Figure 4a presents a comparison 
of a standard GqWo calculation of the M 2 structure us¬ 
ing a grid of 30 frequency points with DFT data. While 
some splitting of the DOS at the Fermi level is evident 
in comparison to DFT, a peak is still observed dX Ep- 
The Go Ho band structure confirms that the while there 


is splitting of the bands, states are evident at the Fermi 
level in the Z — V — A directions and at L. This splitting 
is suggestive of the lower Hubbard band-quasiparticle 
peak-upper Hubbard band splitting^ observed in DMFT 
studies of correlated metals such as paramagnetic V 2 03 r 
however closer inspection reveals a more practical way to 
regard these features. 

Given that the structure contains a Peierls chain, it 
is expected that there will be some bonding-antibonding 
splitting observed, as per the Mi structure. However, as 
this chain does not undergo antiferroelectric distortion, 
the vanadium and oxygen orbitals along the 2 :-axis of the 
chain are not Peierls paired. Therefore a non-bonding, 
metallic band is expected to exist at Ep- Figures 4b-d il¬ 
lustrate this with charge density isosurfaces of the lower 
valence band (Figure 4b), the quasiparticle band (Fig¬ 
ure 4c), and the conduction band (Figure 4d). Clearly, 
the valence band contains all of the bonding density, 
while the quasiparticle and conduction bands are non¬ 
bonding/antibonding. Projecting the density of states 
onto atomic-like orbitals on the Peierls chain vanadium 
atoms (Figure 4e) illustrates that the quasiparticle peak 
is indeed mostly non-bonding 3(1^2 states. The anti¬ 
ferromagnetic (AF) chain (Figure 4f) in contrast exhibits 
more mixed character at Ep. Thus, Figure 4 indicates 
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that in contrast to the Mi structure, standard Gq^o cal¬ 
culations predict that the M 2 structure is metallic due 
to the reduced bonding/antibonding splitting brought 
about by the change in structure. This is at odds with 
the experimentally determined insulating behavior. 

Static PS-G 4 IT 4 calculations however suggest far more 
localised behaviour. The density of states dX Ep (Fig¬ 
ure 5a) is considerably reduced, with the spectral weight 
splitting in a fashion similar to the CuO data of Fig 2 c, 
with some states moving down to form a broad lower 
Hubbard band with the Peierls bonding states, while 
some move upwards into the antibonding band, creat¬ 
ing an upper Hubbard band, separated from the lower 
Hubbard band by a gap of ^ 1 eV. The band structure 
(Figure 5a) reveals considerable depletion of the states 
at the Fermi level, with the valence states shifting down¬ 
wards as the gap opens, and the oxygen p-states at ^ -3 
to -4 eV shift concurrently. The projected DOSs for the 
Peierls chain (Figure 5b) illustrates that the non-bonding 
3 (i^ 2_^2 states observed in the GqWq calculation have 
disappeared from the gap, as have the mixed d-states 
of the AF chain (Figure 5d). It is evident that increas¬ 
ing the effect of electron correlations splits the metallic 
non-bonding states into Hubbard bands, as expected of 
a Mott insulator. 


IV. CONCLUSION 

In summary, accounting for scattering due to on-site 
repulsion in the low q limit of the polarizability in GW 
calculations of Mott insulators correctly reproduces some 
of the signatures of Mott physics. Specifically, for low- 
symmetry structures in which significant Fermi surface 
nesting is not expected, the approach generates the split¬ 
ting of the states at the Fermi level into the upper and 
lower Hubbard bands, and the giant transfer of spec¬ 
tral weight with electron doping. Upon application of 
this extension to the M 2 form of VO 2 it is evident that 
the Ml and M 2 forms are fundamentally different types 
of insulator. The Mi form opens a gap from bonding 
anti-bonding splitting; the bonding valence band wave- 
functions are stabilized with respect to the valence band 
states. This reduces polarizability and unscreens the con¬ 
duction band, giving rise to a correlated band insulator 
This type of insulator is well described by conventional 
Go Wo calculations. The M 2 structure however opens a 
gap due to strong local k-space correlations, it is a Mott 
insulator 

The approach described however rests upon the as¬ 
sumption that after construction of the input wavefunc- 
tions and eigenvalues by DFT, the low q dielectric re¬ 
sponse contains all of the polarizations that create dou¬ 
ble occupancies, and nothing else. While mathematically 
there is certainly a low q dependence due to the pho¬ 
ton propagator, this is not a rigorous enough, or general 
enough approach to be widely applicable. However, the 
computation of the response function is far less taxing 


a 



b 



FIG. 6 . a) Effect of self-consistency iterations on the spin- 
resolved PS-GW data for CuO, b) effect of self-consistency 
on the spin-independent and PS-GW data for CuO (the lower 
computational requirements allow more iterations to be used), 
the data for 4 , 5 and 6 iterations are all plotted in black as 
they practically overwrite each other, and c) effect of self- 
consistency iterations for the M2 VO2 structure. 


than the self-energy, and the possibility exists to imple¬ 
ment a more rigorous approach to the determination of 
transitions which generate double occupancies. If such 
an approach could be found, modified self-energy calcula¬ 
tions of the kind presented here could provide significant 
insight into the nature of Mott insulating materials from 
an ab initio perspective, and thus significantly facilitate 
materials design. 


V. APPENDIX 

Figures 6 a-c present convergence tests using the to¬ 
tal Density of States for GuO spin-resolved (Figure 6 a) 
and spin-independent (Figure 6 b) calculations, and spin- 
independent M 2 VO 2 calculations (Figure 6 c). In all 
cases convergence for static calculations {i.e. those used 
in this study) occurred quickly. For CuO, three itera¬ 
tions was sufficient, while for M 2 VO 2 convergence was 
reached within four iterations. 
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